home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 10 - 1994 / 10.11 Nov 94 / Programmers' Challenge / Factor64.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-06-13  |  21.0 KB  |  1,138 lines  |  [TEXT/KAHL]
  1. /*------------------------------------------------------*/
  2. /* Factor the product N of two primes each of ≤ 32 bits    */
  3. /*------------------------------------------------------*/
  4.  
  5. #include "Factor64.h"
  6.  
  7. /*------------------------------------------------------*/
  8. /*    Factor64() is a simplified version of the quadratic */
  9. /*     sieve using multiple polynomials                    */
  10. /*------------------------------------------------------*/
  11.  
  12. void Factor64(ulong nL,ulong nH,ulong *p1,ulong *p2) {
  13.  
  14.     register ushort ls,k;
  15.     register ulong  i,j;
  16.     
  17.     ushort    p,s,qi,hi,r,c,sP,sN;
  18.     ushort  ts,b,m,bm,br,m2,S[5];
  19.     ulong    d,e,sX,sQ,sq;
  20.  
  21.     Int         B,C,Q,R,T,X,Y;
  22.     ushort  *Ptr,*pf,*sf,*lg,*r1,*r2;
  23.     ushort     *lv,**hm,*hp,**gm,*gp,*gi,*pv,**Xv,*Yv;
  24.  
  25. /*------------------------------------------------------*/
  26. /*    Check whether N is square                            */
  27. /*------------------------------------------------------*/
  28.     
  29.     N[1] = nL & 0xffff; N[2] = nL >> 16;
  30.     N[3] = nH & 0xffff; N[4] = nH >> 16;
  31.     k = 4; while (N[k] == 0) k--; N[0] = k;
  32.  
  33.      sq = floor(sqrt(nL + 4294967296.0*nH));
  34.     
  35.     Set(S,sq); Mul(S,S,S);
  36.  
  37.     if (Comp(S,N) == 0) {
  38.         *p1 = sq;
  39.         *p2 = sq;
  40.         return;
  41.     }
  42.      
  43. /*------------------------------------------------------*/
  44. /*    Check N for small factors (up to 0x800)                */
  45. /*------------------------------------------------------*/
  46.  
  47.     for (k = 0; k < 616; k += 2) {
  48.         p = Prm[k];
  49.         if (Modq(N,p) == 0) { 
  50.             Divq(N,p);
  51.             *p1 = p;
  52.             *p2 = Unset(N);
  53.             return;
  54.         }
  55.     } 
  56.     
  57. /*------------------------------------------------------*/
  58. /*  Allocate memory                                        */
  59. /*------------------------------------------------------*/    
  60.     
  61.     Ptr = malloc(0x20000);
  62.     if (Ptr == 0) {
  63.         printf(" malloc failed \n");
  64.         exit(0);
  65.     }
  66.     X     = (ushort *) Ptr;        Ptr += 20;
  67.     Y     = (ushort *) Ptr;        Ptr += 20;
  68.     B     = (ushort *) Ptr;        Ptr += 20;
  69.     C     = (ushort *) Ptr;        Ptr += 20;
  70.     Q     = (ushort *) Ptr;        Ptr += 20;
  71.     R     = (ushort *) Ptr;        Ptr += 20;
  72.     T     = (ushort *) Ptr;        Ptr += 20;
  73.     
  74.     pf       = (ushort *) Ptr;       Ptr += 120;
  75.     sf       = (ushort *) Ptr;       Ptr += 120;
  76.     lg       = (ushort *) Ptr;       Ptr += 120;
  77.     r1       = (ushort *) Ptr;       Ptr += 120;
  78.     r2       = (ushort *) Ptr;       Ptr += 120;    
  79.     lv       = (ushort *) Ptr;       Ptr += 12000;
  80.  
  81.     hm    = (ushort **) Ptr;    Ptr += 240;
  82.     hm[0] = (ushort *)  Ptr;    Ptr += 14400;        
  83.     gm    = (ushort **) Ptr;    Ptr += 240;
  84.     gm[0] = (ushort *)  Ptr;    Ptr += 24000;
  85.     
  86.     pv       = (ushort *)  Ptr;       Ptr += 120;
  87.     Yv    = (ushort *)  Ptr;    Ptr += 120;
  88.     Xv    = (ushort **) Ptr;    Ptr += 240;
  89.       Xv[0] = (ushort *)  Ptr;                         
  90.             
  91.     for (k = 1; k < 120; k++) {
  92.         hm[k] = hm[k-1] + 120;
  93.          gm[k] = gm[k-1] + 120 + k;
  94.          Xv[k] = Xv[k-1] + 6;
  95.      }
  96.  
  97. /*------------------------------------------------------*/
  98. /*    Parameters for the sieve : ts and qs are cutoffs    */
  99. /*    b is size of prime base and m2 is size of sieve        */
  100. /*------------------------------------------------------*/
  101.  
  102.     k = Bitsize(N);
  103.     
  104.     b = 2 + (5*k)/4;
  105.     bm = b + 3;
  106.     
  107.     if (k > 40) m =  50*k +  400;
  108.     else         m = 100*k - 1600;
  109.      m2 = m + m; 
  110.     
  111.     if (k > 40) ts = 196 + k/2;
  112.     else        ts = 11*k - 24 - (k*k)/8;
  113.     
  114.      sq = ceil(sqrt(sq/m));
  115.      
  116. /*------------------------------------------------------*/
  117. /*    Construct the prime base                            */
  118. /*------------------------------------------------------*/
  119.  
  120. L0:    k = 2;
  121.     i = 0;
  122.     while (k < b) {
  123.         p = Prm[i];
  124.         i++;
  125.         s = Modq(N,p);
  126.         if (qrs(s,p) == 1) {
  127.             pf[k] = p;
  128.             sf[k] = mrt(s,p);
  129.             lg[k] = Prm[i];
  130.             k++;
  131.         }
  132.         i++;
  133.     }
  134.     pf[1] = 2;
  135.  
  136. /*------------------------------------------------------*/
  137. /*    Construct quadratic polynomials as needed            */
  138. /*------------------------------------------------------*/
  139.     
  140.     while (Prm[i] < sq) i += 2;
  141.     hi = i;
  142.     qi = 0;    
  143.  
  144. L1: do {if (hi < 616) sP = Prm[hi];
  145.            else           sP = npr(sP);
  146.            while ((sP&3) == 1) {
  147.              hi += 2;
  148.             if (hi < 616) sP = Prm[hi];
  149.                else           sP = npr(sP);
  150.         }
  151.         sN = Modq(N,sP);
  152.         hi += 2;
  153.     } while (qrs(sN,sP) == -1);
  154.  
  155.     Set(T,sN);
  156.      Dif(T,T,N); 
  157.      Divq(T,sP); 
  158.      d = Modq(T,sP); 
  159.      d *= sP;
  160.      d += sN;
  161.     sX = hrt(d,sP);
  162.     Set(X,sX);
  163.  
  164.     e = (ulong) sP*sP;
  165.     Set(B,e);
  166.     Mulq(B,m);  Dif(B,B,X);
  167.     Mul(C,B,B); Dif(C,C,N);
  168.     Divq(C,sP);    Divq(C,sP);
  169.  
  170. /*------------------------------------------------------*/
  171. /*    Do the sieve for current polynomial                    */
  172. /*------------------------------------------------------*/
  173.  
  174.       if ((B[1] & 1) == 0) {i = 1; j = 0;}
  175.       else                  {i = 0; j = 1;}
  176.  
  177.       if      ((N[1] & 7) == 0) ls = 21;
  178.       else if ((N[1] & 3) == 0) ls = 14;
  179.            else                    ls =  7;
  180.  
  181.      while (i < m2) {
  182.          lv[i] = ls; i += 2;
  183.         lv[j] =  0; j += 2;
  184.     }
  185.  
  186.     for (k = 2; k < b; k++)    {
  187.         p = pf[k]; 
  188.         s = sf[k]; 
  189.         e = sX % p; 
  190.         d = sP % p; 
  191.         d *= d;
  192.         d %= p;
  193.         d = inv(d,p);
  194.  
  195.         if (e < s) e += p;
  196.          i = e-s;
  197.         i *= d;
  198.         i %= p;
  199.         i = m-i;
  200.         i %= p;
  201.         r1[k] = i;
  202.          
  203.         j = e+s;
  204.         j *= d;
  205.         j %= p;
  206.         j = m-j;
  207.         j %= p;
  208.         r2[k] = j;
  209.         
  210.         if (i > j) {
  211.             e = i; 
  212.             i = j; 
  213.             j = e;
  214.         }
  215.  
  216.            ls = lg[k];
  217.           while (j < m2) {
  218.               lv[i] += ls; i += p;
  219.               lv[j] += ls; j += p;
  220.           }
  221.           if (i < m2) lv[i] += ls;
  222.     }
  223.    
  224. /*------------------------------------------------------*/
  225. /*    Factor polynomial to find rows of matrix hm[qi,k]    */
  226. /*------------------------------------------------------*/
  227.  
  228.     for (i = 0; i < m2; i++) if (lv[i] > ts) {
  229.           hp = hm[qi];
  230.           d = (ulong) i*sP;
  231.         Set(T,d);
  232.         Mul(Q,T,T);  Add(Q,Q,C);
  233.         R[0] = B[0]; R[1] = B[1];
  234.         R[2] = B[2]; R[3] = B[3];
  235.         s = i+i; Mulq(R,s);
  236.         if (Comp(Q,R) == 1) hp[0] = 0;
  237.         else                 hp[0] = 1;
  238.           Dif(Q,Q,R);
  239.  
  240.         hp[1] = 0;
  241.          while ((Q[1] & 1) == 0) {
  242.             hp[1] += 1;
  243.                Shiftr(Q);
  244.           }
  245.  
  246.         k = b;
  247.         while (Q[0] > 2) {
  248.             k--;
  249.             if (k == 1) goto L2;
  250.               p = pf[k]; 
  251.                j = i % p;
  252.               if (j == r1[k] || j == r2[k]) {
  253.                   hp[k] = 1;
  254.                   Divq(Q,p);
  255.                   while (Modq(Q,p) == 0) {
  256.                       hp[k] += 1; 
  257.                       Divq(Q,p);
  258.                   }
  259.               }
  260.                else hp[k] = 0;
  261.          }
  262.           sQ = Unset(Q);
  263.           while (sQ > 1) {
  264.               k--;
  265.               if (k == 1) goto L2;
  266.               p = pf[k]; 
  267.                j = i % p;
  268.               if (j == r1[k] || j == r2[k]) {
  269.                   hp[k] = 1;
  270.                   sQ /= p;
  271.                   while (sQ%p == 0) {
  272.                        hp[k] += 1;
  273.                        sQ /= p;
  274.                  }
  275.               }
  276.              else hp[k] = 0;
  277.          }
  278.  
  279.           while (k > 2) {
  280.               k--; 
  281.               hp[k] = 0;
  282.           }
  283.           Mulq(T,sP); 
  284.           Dif(Xv[qi],T,B);
  285.           Yv[qi] = sP;
  286.           qi += 1;
  287.            if (qi == bm) goto L3;
  288. L2:        ;
  289.       }
  290.     goto L1;
  291.  
  292. /*------------------------------------------------------*/
  293. /*    Row reduce gm = hm % 2 and find X^2 = Y^2 (mod N)    */
  294. /*------------------------------------------------------*/
  295.  
  296. L3:    k = N[0];
  297.     g = (double) N[k];
  298.     if (d > 1) g += N[k-1]/65536.0;
  299.     if (d > 2) g += N[k-2]/4294967296.0;
  300.     if (d > 3) g += 1/4294967296.0;
  301.     
  302.     for (r = 0; r < bm; r++) { 
  303.         hp = hm[r]; gp = gm[r];
  304.         for (c = 0; c < b; c++) 
  305.             gp[c]  = hp[c] & 1;
  306.         br = b + r;
  307.         for (c = b; c < br; c++) gp[c] = 0;    
  308.         gp[br] = 1;
  309.     }
  310.   
  311.     for (r = 0; r < bm; r++) {
  312.         br = b + r;
  313.         gp = gm[r];
  314.  
  315.         c = b - 1;  
  316.         while (gp[c] == 0 && c > 0) c--;
  317.  
  318.         if (c > 0 || gp[0] == 1) { 
  319.             for (i = r+1; i < bm; i++) {
  320.                    gi = gm[i];
  321.                 if (gi[c] == 1) {
  322.                     gi[c] = 0;
  323.                     for (k = 0; k < c ; k++)
  324.                          gi[k] ^= gp[k]; 
  325.                     for (k = b; k <= br; k++)
  326.                          gi[k] ^= gp[k];
  327.                   }
  328.             }
  329.         }
  330.         else {
  331.             for (i = 1; i < b; i++) pv[i] = 0;
  332.  
  333.             Set(X,1); Set(Y,1);
  334.             for (k = b; k <= br; k++) 
  335.                 if (gp[k] == 1) {
  336.                     j = k - b;
  337.                     Mul(X,X,Xv[j]); 
  338.                     if (X[0] > 12) Mod(X);    
  339.                     Mulq(Y,Yv[j]);  
  340.                     if (Y[0] > 12) Mod(Y);    
  341.                     for (i = 1; i < b; i++)
  342.                         pv[i] += (long) hm[j][i];
  343.                 }
  344.                Mod(X);
  345.                
  346.                k = pv[1] >> 1;
  347.                while (k > 15) {
  348.                    for (i = Y[0]; i > 0; i--)
  349.                        Y[i+1] = Y[i];
  350.                    Y[1] = 0; 
  351.                    Y[0] += 1;
  352.                    k -= 16;
  353.                    if (Y[0] > 12) Mod(Y);
  354.                }
  355.                Mulq(Y,1 << k);
  356.  
  357.             for (i = 2; i < b; i++) {
  358.                 j = pv[i];
  359.                 if (j == 0) continue;
  360.                 if (j == 2) Mulq(Y,pf[i]);
  361.                 else {
  362.                     T[1] = pf[i]; 
  363.                     T[0] = 1;
  364.                     while (j > 2) {
  365.                         j -= 2; 
  366.                         Mulq(T,pf[i]);
  367.                     }
  368.                     Mul(Y,Y,T);    
  369.                 }
  370.                 if (Y[0] > 12) Mod(Y);
  371.             }
  372.             Mod(Y);
  373.             
  374.             Dif(T,X,Y);    
  375.             Add(R,X,Y);
  376.             if (Comp(N,R) <= 0) Dif(R,R,N);
  377.             if (T[0] == 0 || R[0] == 0) continue;
  378.             *p1 = Gcd(T); 
  379.             *p2 = Gcd(R);
  380.             return;    
  381.           }
  382.     } 
  383.     bm += 5; 
  384.     ts -= 2;
  385.     goto L0;     
  386. }
  387.  
  388. /*------------------------------------------------------*/
  389. /*  End of the function Factor64().                     */
  390. /*------------------------------------------------------*/
  391.  
  392. /*------------------------------------------------------*/
  393. /*    Inverse of b modulo m  (Euclid’s algorithm)            */
  394. /*------------------------------------------------------*/
  395.  
  396. ulong inv(ulong b,ushort m) {
  397.  
  398.     register ulong u,v,t,n;                                
  399.  
  400.     u = 1;
  401.     v = 0;
  402.     t = 1;
  403.     n = m;
  404.     
  405.     for (;;) {        
  406.         while ((b & 1) == 0) {            
  407.             v <<= 1;                    
  408.             if (t & 1) t += m;                
  409.             t >>= 1;                    
  410.             b >>= 1;
  411.         }
  412.         if (b == 1) break;                        
  413.     
  414.         if (b > n) {
  415.             b -= n;    
  416.             u += v;                        
  417.         }
  418.         else {
  419.             n -= b;                
  420.             v += u;                    
  421.         }
  422.     
  423.         while ((n & 1) == 0) {
  424.             u <<= 1;                        
  425.             if (t & 1) t += m;            
  426.             t >>= 1;                    
  427.             n >>= 1;
  428.         }
  429.     }
  430.     t *= u;
  431.     t %= m;
  432.          
  433.     return t;
  434. }
  435.  
  436. /*------------------------------------------------------*/
  437. /*    Determine if n is square modulo p  (QR algorithm)   */
  438. /*------------------------------------------------------*/
  439.  
  440. short qrs(ushort n,ushort p) {
  441.         
  442.      register short  j;
  443.      register ushort n2,n4;
  444.  
  445.      if (n == 1) return 1;
  446.  
  447.      j = 1;
  448.      for (;;) {
  449.           n2 = n + n;
  450.           n4 = n2 + n2;
  451.          p %= n4;
  452.          while (p > n) {
  453.                if (n & 2)  j = -j;
  454.                if (p > n2) p -= n2;
  455.                else        p = n2 - p;
  456.          }
  457.          if (p == 1) return j;
  458.          n %= p;
  459.          if (n == 1) return j;
  460.       }
  461. }
  462.  
  463. /*------------------------------------------------------*/
  464. /*    Square root of n modulo p                            */
  465. /*------------------------------------------------------*/
  466.  
  467. ushort mrt(ushort n,ushort p) {
  468.  
  469.        register ulong q,s;
  470.        register long  x,u,v;
  471.        ulong m;
  472.       long  t,r;
  473.  
  474.     if (n == 1) return 1;
  475.     q = (p+1) >> 1;
  476.  
  477.      m = n;
  478.      if ((q & 1) == 0) {
  479.         q >>= 1;
  480.         s   = 1;
  481.         while (q > 0) {
  482.             if (q & 1) {
  483.                 s *= m;
  484.                 s %= p;
  485.             }
  486.             m *= m;
  487.             m %= p;
  488.             q >>= 1;
  489.         }
  490.            if (s+s > p) s = p-s;
  491.          return s;
  492.        }
  493.  
  494.       s = 0;
  495.     t = n;
  496.     while (qrs(t,p) == 1) {
  497.            s += 1;
  498.            t += 1-s-s;
  499.           if (t%p == 0) {
  500.               if (s+s > p) s = p-s;
  501.              return s;
  502.          }
  503.          if (t < p) t += p;
  504.       }
  505.  
  506.     q >>= 1;
  507.     n = t;
  508.     r = s;
  509.     u = 1;
  510.     v = 1;
  511.     while (q > 0) {
  512.         x = n*u;
  513.         x %= p;
  514.         x *= u;
  515.         x %= p;
  516.         t = r*r;
  517.         t %= p;
  518.         if (t >= x) x = t-x;
  519.         else        x = t-x+p;
  520.         u *= r;
  521.         u %= p;
  522.         u += u;
  523.         u %= p;
  524.         r = x;
  525.            if (q & 1) {
  526.               x = n*u;
  527.               x %= p;
  528.             x *= v;
  529.             x %= p;
  530.             t = s*r;
  531.             t %= p;
  532.             if (t >= x) x = t-x;
  533.             else        x = t-x+p;
  534.             v *= r;
  535.             v %= p;
  536.             v += s*u;
  537.             v %= p;
  538.             s = x;
  539.            }
  540.            q >>= 1;
  541.       }
  542.  
  543.     if (s+s > p) s = p-s;
  544.        return s;
  545. }
  546. /*------------------------------------------------------*/
  547. /*    Square root of n modulo p^2 for p = 3 (mod 4)         */
  548. /*------------------------------------------------------*/
  549.  
  550. ulong hrt(ulong n,ushort p) {
  551.  
  552.     register ulong  s,t,x,y;
  553.     ulong  m;
  554.  
  555.     m = n%p;
  556.     s = m;
  557.     y = 1;
  558.     t = p-3;
  559.     t >>= 2;
  560.     while (t > 0) {
  561.           if (t&1) {
  562.                y *= s;
  563.                y %= p;
  564.            }
  565.           s *= s;
  566.           s %= p;
  567.            t >>= 1;
  568.       }
  569.     x = m*y;
  570.     x %= p;
  571.  
  572.     s = (ulong) p*p;
  573.     t = x*x;
  574.     if (n < t) n += s;
  575.     n -= t;
  576.     n /= p;
  577.     n *= y;
  578.     n %= p;
  579.  
  580.     t = p;
  581.     t += 1;
  582.     t >>= 1;
  583.     n *= t;
  584.     n %= p;
  585.  
  586.     n *= p;
  587.     n += x;
  588.     if (n+n > s) n = s-n;
  589.     return n;
  590. }
  591.  
  592. /*------------------------------------------------------*/
  593. /*    Get next prime                                          */
  594. /*------------------------------------------------------*/
  595.  
  596. ushort npr(ushort p) {
  597.  
  598.     register ushort d,s,k;
  599.  
  600.       do {p += 2;
  601.           s = floor(sqrt(p));
  602.           k = 0;
  603.           d = 3;
  604.         while (d <= s) {
  605.              if (p%d == 0) break;
  606.             k += 2;
  607.             d = Prm[k];
  608.         }
  609.      }
  610.      while (d <= s);
  611.      return p;
  612. }
  613.  
  614. /*------------------------------------------------------*/
  615. /*    Addition  S = A + B                                    */
  616. /*------------------------------------------------------*/
  617.  
  618. void Add(Int S,Int A,Int B) {
  619.  
  620.     register ushort *pH, *pL;
  621.     register ulong  t;
  622.     register ushort c,k;
  623.     ushort s,dH,dL;
  624.     
  625.     if (A[0] > B[0] ) {
  626.         pH = A; dH = A[0];
  627.         pL = B; dL = B[0];
  628.     }
  629.     else {
  630.         pH = B; dH = B[0];
  631.         pL = A; dL = A[0];
  632.     }
  633.     if (dL == 0) {
  634.         if (S != pH) 
  635.             for (k = 0; k <= dH; k++) S[k] = pH[k];
  636.         return;
  637.     }
  638.     
  639.     k = 0;
  640.     c = 0;
  641.     while (k < dL) { 
  642.         k++;
  643.         t = (ulong) pH[k] + pL[k] + c;
  644.         if (t >= 0x10000) {
  645.             t -= 0x10000; 
  646.             c = 1;
  647.         }
  648.         else c = 0;
  649.         S[k] = t;
  650.     }
  651.     while (c == 1 && k < dH) {
  652.         k++;
  653.         s = pH[k];
  654.         if (s == 0xFFFF) {
  655.             S[k] = 0; 
  656.             c = 1;
  657.         }
  658.         else {
  659.             S[k] = s + 1;
  660.             c = 0;
  661.         }
  662.     }
  663.     while (k < dH) {
  664.         k++;
  665.         S[k] = pH[k];
  666.     }
  667.     if (c == 1) {
  668.         dH += 1;
  669.         S[dH] = 1;
  670.     }
  671.     S[0] = dH;                    
  672. }
  673.  
  674. /*------------------------------------------------------*/
  675. /*    Difference  D = |A - B|                                */
  676. /*------------------------------------------------------*/
  677.  
  678. void Dif(Int D,Int A,Int B) {
  679.  
  680.     register ushort *pH, *pL;
  681.     register long   t;
  682.     register ushort c,k;
  683.     ushort   s,dH,dL;
  684.     short    e;
  685.  
  686.     k = A[0];
  687.     if (k > B[0]) e = 1;
  688.     else {
  689.         if (k < B[0]) e = -1;
  690.         else {
  691.             while (A[k] == B[k] && k > 0) k--;
  692.             if (k == 0) {
  693.                 D[0] = 0;
  694.                 return;
  695.             }
  696.             if (A[k] > B[k]) e = 1;
  697.             else e = -1;
  698.         }
  699.     }
  700.     if (e == 1) {
  701.         pH = A; dH = A[0];
  702.         pL = B; dL = B[0];
  703.     }
  704.     else {
  705.         pH = B; dH = B[0];
  706.         pL = A; dL = A[0];
  707.     }
  708.     if (dL == 0) {
  709.         if (D != pH) 
  710.             for (k = 0; k <= dH; k++) D[k] = pH[k];
  711.         return;
  712.     }
  713.  
  714.     c = 0;
  715.     k = 0;
  716.     while (k < dL) {
  717.         k++;        
  718.         t = (long) pH[k] - pL[k] - c;
  719.         if (t < 0) {
  720.             t += 0x10000; 
  721.             c = 1;
  722.         }
  723.         else c = 0;
  724.         D[k] = t;
  725.     }
  726.     while (c == 1 && k < dH) {
  727.         k++;
  728.         s = pH[k];
  729.         if (s == 0) {
  730.             D[k] = 0xFFFF; 
  731.             c = 1;
  732.         }
  733.         else {
  734.             D[k] = s - 1;
  735.             c = 0;
  736.         }
  737.     }
  738.     while (k < dH) {
  739.         k++;
  740.         D[k] = pH[k];
  741.     }
  742.     k = dH;
  743.     while (D[k] == 0 && k > 0) k--;
  744.     D[0] = k;
  745. }
  746.  
  747. /*------------------------------------------------------*/
  748. /*    Multiply  P = A * B                                    */
  749. /*------------------------------------------------------*/
  750.  
  751. void Mul(Int P,Int A,Int B) {
  752.  
  753.     register ushort *pB;
  754.     register ulong  t;
  755.     register ushort s,n,k;
  756.     ushort d,j;
  757.     
  758.     if (A[0] == 0 || B[0] == 0) {
  759.         P[0] = 0;
  760.         return;
  761.     }
  762.     d = A[0] + B[0];
  763.     
  764.     for (k = 1; k <= d; k++) bufm[k] = 0;
  765.     
  766.     pB = B;
  767.     j = 0;
  768.     while (j < A[0]) {
  769.         j++;
  770.         s = A[j];
  771.         k = 0;
  772.         n = j;
  773.         while (k < pB[0]) {
  774.             k++;
  775.             t = (ulong) s * pB[k];
  776.             bufm[n] += t & 0xFFFF;
  777.             n++;
  778.             bufm[n] += t >> 16;
  779.         }        
  780.     }
  781.     
  782.     k = 1;
  783.     while (k < d) {
  784.         t = bufm[k];
  785.         bufm[k] = t & 0xFFFF;
  786.         k++;
  787.         bufm[k] += t >> 16;
  788.     }
  789.     
  790.     if (bufm[d] == 0) d -= 1;    
  791.     
  792.     for (k = 1; k <= d; k++) P[k] = bufm[k];
  793.     P[0] = d;
  794. }
  795.  
  796. /*------------------------------------------------------*/
  797. /*    Quick multiply  A = A * b                           */
  798. /*------------------------------------------------------*/
  799.  
  800. void Mulq(Int A,ushort b) {
  801.  
  802.     register ushort *pA;
  803.     register ulong  t,w,c;
  804.     register ushort d,k;
  805.     
  806.     d = A[0];
  807.     if (d == 0) return;
  808.     
  809.     pA = A;
  810.     if (b == 0) {
  811.         pA[0] = 0;
  812.         return;
  813.     }
  814.     
  815.     c = 0;
  816.     k = 0;
  817.     while (k < d) {
  818.         k++;
  819.         t = (ulong) pA[k] * b;
  820.         w = (t & 0xFFFF) + c;
  821.         c = t >> 16;
  822.         if (w >= 0x10000) {
  823.             w -= 0x10000;
  824.             c += 1;
  825.         }
  826.         pA[k] = w;
  827.     }
  828.     if (c > 0) {
  829.         d += 1;
  830.         pA[d] = c;
  831.     }
  832.     A[0] = d;
  833. }
  834.  
  835. /*------------------------------------------------------*/
  836. /*    Modulo  R = R % N                                    */
  837. /*------------------------------------------------------*/
  838.  
  839. void Mod(Int R) {
  840.  
  841.     register ulong  t;
  842.     register long   w;
  843.     register ushort k,j;
  844.     ushort d,n,tL,tH;
  845.     ulong  c;    
  846.     short  e;
  847.     double f;
  848.  
  849.     d = N[0];
  850.     n = R[0];
  851.  
  852.     if (d > n) return;
  853.  
  854.     for (k = 1; k <= n; k++) buf[k] = R[k];
  855.  
  856.     while (n >= d) {
  857.     
  858.         f = buf[n]*65536.0;
  859.         if (n > 1) f += buf[n-1];
  860.         t = f/g;
  861.     
  862.         e = n - d;
  863.         if (e > 0) {
  864.             tL = t & 0xFFFF;
  865.             tH = t >> 16;
  866.         }
  867.          else {
  868.              tL = t >> 16;
  869.              if (tL == 0) {
  870.                  if (t < 0xFFFF) break;
  871.                  else {
  872.                      k = d;
  873.                      while (buf[k] == N[k] && k > 0) k--;
  874.                      if (k > 0 && buf[k] < N[k]) break;
  875.                      tL = 1;
  876.                  }
  877.              }
  878.              tH = 0;
  879.              e = 1;
  880.          }
  881.  
  882.         c = 0;
  883.         j = e;
  884.         for (k = 1; k <= d; k++) {
  885.             t = (ulong) N[k]*tL + c;
  886.             w = (ulong) buf[j] - (t & 0xFFFF);
  887.             t >>= 16;
  888.             if (w < 0) {
  889.                 buf[j] = 0x10000 + w;
  890.                 t += 1;
  891.             }
  892.             else buf[j] = w;
  893.             j++;
  894.             w = (ulong) buf[j] - t;
  895.             if (w < 0) {
  896.                 buf[j] = 0x10000 + w;
  897.                  c = 0x10000;
  898.             }
  899.             else {
  900.                 buf[j] = w;
  901.                 c = 0;
  902.             }
  903.         }
  904.         
  905.         if (tH > 0) {
  906.             c = 0;
  907.             j = e + 1;
  908.             for (k = 1; k <= d; k++) {
  909.                 t = (ulong) N[k]*tH + c;
  910.                 w = (ulong) buf[j] - (t & 0xFFFF);
  911.                 t >>= 16;
  912.                 if (w < 0) {
  913.                     buf[j] = 0x10000 + w;
  914.                     t += 1;
  915.                 }
  916.                 else buf[j] = w;
  917.  
  918.                 if (k == d) break;    
  919.                 j++;
  920.                 w = (ulong) buf[j] - t;
  921.                 if (w < 0) {
  922.                     buf[j] = 0x10000 + w;
  923.                      c = 0x10000;
  924.                 }
  925.                 else {
  926.                     buf[j] = w;
  927.                     c = 0;
  928.                 }
  929.             }
  930.         }
  931.          while (buf[n] == 0 && n > 0) n--;
  932.          if (n == d && buf[d] < N[d]) break;
  933.     }
  934.  
  935.     for (k = 1; k <= n; k++) R[k] = buf[k];
  936.     R[0] = n;
  937. }
  938.  
  939. /*------------------------------------------------------*/
  940. /*    Quick modulo  A = A % m                                */
  941. /*------------------------------------------------------*/
  942.  
  943. ushort Modq(Int A,ushort m) {
  944.  
  945.     register ulong  z,t;
  946.     ushort n,e;
  947.  
  948.     n = A[0];
  949.     if (n == 0) return 0;
  950.  
  951.     for (e = 1; e <= n;    e++) buf[e] = A[e];
  952.  
  953.     while (n > 1) {
  954.         e = n - 1;
  955.         t = ((ulong) buf[n] << 16) + buf[e];
  956.         z = t/m; 
  957.  
  958.         t -= z*m;
  959.         buf[e] = t;
  960.         if (t == 0) do e--;
  961.         while (buf[e] == 0 && e > 0);
  962.         n = e;
  963.     }
  964.     return buf[1] % m;
  965. }
  966.  
  967. /*------------------------------------------------------*/
  968. /*    Quick divide  A = A / m                                */
  969. /*------------------------------------------------------*/
  970.  
  971. void Divq(Int A,ushort m) {
  972.  
  973.     register ulong  z,t;
  974.     ushort n,e;
  975.  
  976.     n = A[0];
  977.     if (n == 0) return;
  978.  
  979.     for (e = 1; e <= n;    e++) {
  980.         buf[e] = A[e];
  981.         A[e] = 0;
  982.     }
  983.  
  984.     while (n > 1) {
  985.         e = n - 1;
  986.         t = ((ulong) buf[n] << 16) + buf[e];
  987.         z = t/m;
  988.  
  989.         if (z < 0x10000) A[e] = z;
  990.         else {
  991.             A[e] = z & 0xFFFF;
  992.             A[n] = z >> 16;
  993.         }
  994.         
  995.         t -= z*m;
  996.         buf[e] = t;
  997.         if (t == 0) do e--;
  998.         while (buf[e] == 0 && e > 0);
  999.         n = e;
  1000.     }
  1001.  
  1002.     n = buf[1];
  1003.     if (n >= m) A[1] = n/m;
  1004.     if (A[A[0]] == 0) A[0] -= 1;
  1005. }
  1006.  
  1007. /*------------------------------------------------------*/
  1008. /*    Shift right  X = X >> 1                                */
  1009. /*------------------------------------------------------*/
  1010.  
  1011. void Shiftr(Int X) {
  1012.  
  1013.     register ulong  t;
  1014.     register ushort i,j,c,d;
  1015.  
  1016.     d = X[0];
  1017.     if (d == 0) return;
  1018.  
  1019.     i = 0;
  1020.     j = 1;
  1021.     c = X[1] >> 1;
  1022.     while (j < d) {
  1023.         j++;
  1024.         t = (ulong) X[j] << 15;
  1025.         t += c;
  1026.         i++;
  1027.         X[i] = t & 0xFFFF;
  1028.         c = t >> 16;
  1029.     }
  1030.     if (c == 0) d -= 1;
  1031.     else X[d] = c;
  1032.     X[0] = d;
  1033. }
  1034.  
  1035. /*------------------------------------------------------*/
  1036. /*    Compare : {+1  0  -1} as {X > Y   X == Y   X < Y}    */
  1037. /*------------------------------------------------------*/
  1038.  
  1039. short Comp(Int X,Int Y) {
  1040.  
  1041.     register ushort d;
  1042.  
  1043.     d = X[0];
  1044.  
  1045.     if (d > Y[0]) return  1;
  1046.     if (d < Y[0]) return -1;
  1047.  
  1048.     while (X[d] == Y[d] && d > 0) d--;
  1049.  
  1050.     if (d == 0) return 0;
  1051.     if (X[d] > Y[d]) return 1;
  1052.     return -1;
  1053. }
  1054.  
  1055. /*------------------------------------------------------*/
  1056. /*    Convert from unsigned long to Int                    */
  1057. /*------------------------------------------------------*/
  1058.  
  1059. void Set(Int X, ulong n) {
  1060.  
  1061.     if (n == 0) {
  1062.         X[0] = 0;
  1063.         return;
  1064.     }
  1065.     if (n < 0x10000) {
  1066.         X[0] = 1;
  1067.         X[1] = n;
  1068.         return;
  1069.     }
  1070.     X[0] = 2;
  1071.     X[1] = n & 0xffff;
  1072.     X[2] = n >> 16;
  1073. }
  1074.  
  1075. /*------------------------------------------------------*/
  1076. /*    Convert from Int to unsigned long                    */
  1077. /*------------------------------------------------------*/
  1078.  
  1079. ulong Unset(Int X) {
  1080.  
  1081.     register ulong  n;
  1082.     register ushort d;
  1083.  
  1084.     d = X[0];
  1085.     if (d == 0) return 0;
  1086.     n = (ulong) X[1];
  1087.     if (d == 1) return n;
  1088.     return (n + ((ulong) X[2] << 16));
  1089. }
  1090.  
  1091. /*------------------------------------------------------*/
  1092. /*    Number of Bits in X                                    */
  1093. /*------------------------------------------------------*/
  1094.  
  1095. ulong Bitsize(Int X) {
  1096.  
  1097.     register ushort d,t;
  1098.     register ulong  n;
  1099.  
  1100.     d = X[0];
  1101.     if (d == 0) return 0;
  1102.  
  1103.     n = (ulong) d << 4;
  1104.     t = 0x8000;
  1105.     while ((t & X[d]) == 0) {
  1106.         n  -= 1;
  1107.         t >>= 1;
  1108.     }
  1109.     return n;
  1110. }
  1111.  
  1112. /*------------------------------------------------------*/
  1113. /*    The greatest common divisor of A and N                */
  1114. /*------------------------------------------------------*/
  1115.  
  1116. ulong Gcd(Int A) {
  1117.  
  1118.     register long k;
  1119.  
  1120.     for (k = 0; k < 5; k++) buf[k] = N[k];
  1121.  
  1122.     while ((A[1]&1) == 0) Shiftr(A);
  1123.  
  1124.     for (;;)
  1125.         switch (Comp(A,buf)) {
  1126.         case  1 :      Dif(A,A,buf);
  1127.                       while ((A[1]&1) == 0) Shiftr(A);
  1128.                       break;
  1129.         case -1 :      Dif(buf,buf,A);
  1130.                       while ((buf[1]&1) == 0) Shiftr(buf);
  1131.                       break;
  1132.         case  0 :     return Unset(A);
  1133.         }
  1134. }
  1135.  
  1136. /*------------------------------------------------------*/
  1137. /*    End of file Factor64.c                                */
  1138. /*------------------------------------------------------*/